#ifndef _FEMSPACE_DEALII_UNSTRUCT_H_
#define _FEMSPACE_DEALII_UNSTRUCT_H_
#include "AMGPCGsolver.h"
#include "FEMSpace.h"
#include <cassert>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
// using DealiiSpMat = dealii::SparseMatrix<double>;
using DealiiRHS = dealii::Vector<double>;
using DealiiSolution = dealii::Vector<double>;
/**
 * @brief
 * 返回map，存放边界上自由度的全局标号，和边界函数(仅Dirichlet)在对应自由度上的函数值。
 *
 * @param _boundary_function 边界函数
 * @param _boundary_dofs 边界上的自由度
 */
template <unsigned int DIM>
std::map<int, double> interpolateBoundaryValues_dealii(
    EasyMesh *easy_mesh, vector<function<double(double *)>> &_boundary_function,
    std::set<Dofs<DIM>> &_boundary_dofs) {

    std::map<int, double> boundary_values;
    for (auto &dofs : _boundary_dofs) {
        double p[2];
        for (int i = 0; i < 2; ++i) {
            p[i] = dofs[i];
        }
        // cout << "自由度信息:" << dofs << endl;
        // cout << "getGlobalindex\n";
        int index = dofs.getGlobalIndex();
        // EasyMesh中根据自由度index取出boundarymark
        // cout << "getBoundaryMark:\n";
        int boundarymark = easy_mesh->getBoundaryMark(index);
        assert(boundarymark > 0);
        // cout <<"get value\n";
        // cout <<"boundary mark is" << boundarymark << endl;
        double value = _boundary_function[boundarymark - 1](p);
        // cout <<"insert into boundary_values\n";
        boundary_values.insert({index, value});
    }
    return boundary_values;
};

/**
 * @brief 对刚度矩阵和右端项应用Dirichlet条件
 *
 * @param system_matrix 刚度矩阵
 * @param system_rhs 右端项
 */
void applyBoundaryValues_dealii(std::map<int, double> &boundary_values,
                                DealiiSpMat &system_matrix,
                                DealiiRHS &system_rhs,
                                const SparsityPattern &sparsity_pattern) {
    // const SparsityPattern &sparsity_pattern =
    //     system_matrix.get_sparsity_pattern();
    for (auto &it : boundary_values) {
        // 边界自由度的全局编号。
        int k = it.first;
        // 对应自由的边界函数值。
        double bnd_value = it.second;
        // 第k行的非零元个数
        int l = sparsity_pattern.row_length(k);
        system_rhs[k] = bnd_value * system_matrix(k, k);
        for (int m = 0; m < l; m++) {
            int col = sparsity_pattern.column_number(k, m);
            if (col == k)
                continue;
            system_rhs[col] -= system_matrix(k, col) * bnd_value;
            system_matrix.set(k, col, 0);
            system_matrix.set(col, k, 0);
        }
    }
}

class Possion_2D_Dealii : public FEMSpace<2> {
  public:
    Possion_2D_Dealii();
    Possion_2D_Dealii(EasyMesh *easy_mesh, TemplateElement<2> *_element,
                      Equation<2> *_equation, Domain<2> *_domain);
    Possion_2D_Dealii(const Possion_2D_Dealii &p);
    ~Possion_2D_Dealii(){/*delete _mesh;delete _element;delete _equation;*/};
    virtual void setup_system();
    virtual void assemble_system();
    virtual void dealwithBoundaryCondition();
    virtual void solve();

    double LocalStiffMatrix(int i, int j);
    const DealiiSolution& getSolution() const;
    const DealiiSpMat& A() const;
    const DealiiRHS& Rhs() const;

  protected:
    DealiiSpMat system_matrix;
    DealiiRHS system_rhs;
    DealiiSolution solution;
    Domain<2> *domain;
    EasyMesh *mesh;
    TemplateElement<2> *element;
    Equation<2> *equation;
    SparsityPattern sparsity_pattern;
};

double Possion_2D_Dealii::LocalStiffMatrix(int i, int j) {
    double a = 0;
    for (int k = 0; k < element->n_GaussPnt(); k++) {
        double xi = element->GaussionPoint(k)[0];
        double eta = element->GaussionPoint(k)[1];
        double b = 0;
        for (int l = 0; l < element->gradient(xi, eta, i).size(); l++)
            b += element->gradient(xi, eta, i)[l] *
                 element->gradient(xi, eta, j)[l];
        a = a + element->det_Jacobi(xi, eta) * element->GaussionWeight(k) * b;
    }
    return a;
}

const DealiiSolution& Possion_2D_Dealii::getSolution() const { return solution; }

Possion_2D_Dealii::Possion_2D_Dealii() {
    // int num = mesh->getTotalNumDofs();
    // system_matrix = DealiiSpMat(num, num);
    // system_rhs = DealiiRHS(num);
}
Possion_2D_Dealii::Possion_2D_Dealii(EasyMesh *easy_mesh,
                                     TemplateElement<2> *_element,
                                     Equation<2> *_equation,
                                     Domain<2> *_domain) {
    domain = _domain;
    mesh = easy_mesh;
    element = _element;
    equation = _equation;
    int num = mesh->getTotalNumDofs();
    system_rhs = DealiiRHS(num);
    solution = DealiiRHS(num);
}

Possion_2D_Dealii::Possion_2D_Dealii(const Possion_2D_Dealii &p) {
    domain = p.domain;
    mesh = p.mesh;
    element = p.element;
    equation = p.equation;
    system_matrix = p.system_matrix;
    system_rhs = p.system_rhs;
    sparsity_pattern.copy_from(p.sparsity_pattern);
}

void Possion_2D_Dealii::setup_system() {
    /**
     * set the sizes of system_matrix, right hand side vector ans the solution
     * vector;
     */
    int num = mesh->getTotalNumDofs();
    std::cout << "Number of degrees of freedom: " << num << std::endl;

    sparsity_pattern.reinit(num, num, num);
    for (int index = 0; index < mesh->getNumGrid(); ++index) {

        // bind the TemplateElement with current cell.
        Grid<2> cell = mesh->getGrid(index);
        element->reinit(cell, *mesh, index);
        const unsigned int dofs_per_cell = element->n_Dofs();

        for (int i = 0; i < dofs_per_cell; i++) {
            for (int j = 0; j < dofs_per_cell; j++) {
                sparsity_pattern.add(element->getGlobalIndex(i),
                                     element->getGlobalIndex(j));
            }
        }
    }

    sparsity_pattern.compress();
    system_matrix.reinit(sparsity_pattern);
    // cout << "成员变量sparsity_pattern:\n";
    // sparsity_pattern.print(std::cout);
    // const SparsityPattern &matrix_sp = system_matrix.get_sparsity_pattern();
    // cout << "matrix_sp:\n";
    // matrix_sp.print(std::cout);
    solution.reinit(num);
    system_rhs.reinit(num);
}

void Possion_2D_Dealii::assemble_system() {
    /**
     * Get the right_hand_side function and boundary function
     */
    auto right_hand_side = equation->RightHandsTermFunction();

    /**
     * for each Grid, we construct a cell_matirx and cell_rhs.
     */
    const unsigned int dofs_per_cell = element->n_Dofs();
    DealiiRHS cell_rhs(dofs_per_cell);

    /*assemble the global matrix and vector cell-by-cell.*/
    for (int index = 0; index < mesh->getNumGrid(); ++index) {
        // bind the TemplateElement with current cell.
        Grid<2> cell = mesh->getGrid(index);
        element->reinit(cell, *mesh, index);

        for (int i = 0; i < dofs_per_cell; i++) {
            for (int j = 0; j < dofs_per_cell; j++) {
                system_matrix.add(element->getGlobalIndex(i),
                                  element->getGlobalIndex(j),
                                  LocalStiffMatrix(i, j));
            }
        }

        cell_rhs = DealiiRHS(dofs_per_cell);
        for (int q_index = 0; q_index < element->n_GaussPnt(); ++q_index) {
            Point<2> pnt = element->GaussionPoint(q_index);
            double xi = pnt[0], eta = pnt[1];
            double x = element->Global_x(xi, eta);
            double y = element->Global_y(xi, eta);
            double p[2] = {x, y};
            for (int i = 0; i < dofs_per_cell; ++i) {
                cell_rhs(i) += (element->phi(xi, eta, i) * right_hand_side(p) *
                                element->det_Jacobi(xi, eta) *
                                element->GaussionWeight(q_index));
            }
        }
        for (int i = 0; i < dofs_per_cell; ++i)
            system_rhs[element->getGlobalIndex(i)] += cell_rhs(i);
    }
}

void Possion_2D_Dealii::dealwithBoundaryCondition() {
    vector<function<double(double *)>> boundary_function =
        equation->BoundaryConditionFunction();
    // cout << "boundary function的size为:" << boundary_function.size() << endl;
    std::set<Dofs<2>> boundary_dofs;
    boundary_dofs = mesh->getBoundaryDofs();
    // for (auto dof = boundary_dofs.begin(); dof != boundary_dofs.end(); ++dof) {
    //     cout << *dof << "\n";
    // }
    // cout << endl;
    // TODO:采用新函数计算边界自由度上的函数值。
    // cout << "interpolateBoundaryValue\n" << endl;
    std::map<int, double> boundary_values;
    boundary_values = interpolateBoundaryValues_dealii(mesh, boundary_function,
                                                       boundary_dofs);
    applyBoundaryValues_dealii(boundary_values, system_matrix, system_rhs,
                               sparsity_pattern);
}

void Possion_2D_Dealii::solve() {
    /// @brief 建立AMGPCG
    /// @param _Mat 求解矩阵
    /// @param _rebuild 是否要进行矩阵重建(刚度阵这一项为true, 质量阵为false)
    /// @param _maxitnum CG迭代最大迭代次数
    /// @param _tol 迭代阈值
    /// @param theta 强连接点阈值
    /// @param _smoothstep 每一层网格做的GS次数
    /* AMGPCG(DealiiDealiiSpMat &_Mat , bool _rebuild = true, int _maxitnum =
     * 200, double _tol = 1e-13, double theta = 0.25, int _smoothstep = 1); */

    /// 使用方式
    // AMGPCG AMGPCGsolver(system_matrix,true,20,1e-8,0.25,2);
    // solution = AMGPCGsolver.solve(system_rhs);

    //CG
    dealii::SolverControl  solver_control(1000, 1e-10);
    dealii::SolverCG<dealii::Vector<double>> solver(solver_control);
    solver.solve(system_matrix, solution, system_rhs, dealii::PreconditionIdentity());

    // We have made one addition, though: since we suppress output from the
    // linear solvers, we have to print the number of iterations by hand.
    std::cout << "   " << solver_control.last_step()
      << " CG iterations needed to obtain convergence." << std::endl;
}

const DealiiSpMat& Possion_2D_Dealii::A() const
{
  return system_matrix;
}
const DealiiRHS& Possion_2D_Dealii::Rhs() const
{
  return system_rhs;
}

#else
#endif
